Surface Matching: example on human scapula
Contents
1.1. Surface Matching: example on human scapula¶
This notebook presents a method to match a point cloud taken with any scanning system to a digital 3D scan in stl format.
1.1.1. Introduction¶
1.1.1.1. Method outlines¶
Surface matching is based on an optimization strategy. The objective of this optimization is to identify the parameters of the reference change allowing to pass from the coordinate system of the STL file to the coordinate system of the scan system.
This change of reference frame is formalized by a passage matrix formed by a rotation vector and a translation vector according to the Rodrigues formalism. There are thus a total of 6 DOF.
The optimization parameters of the cost function will therefore be these two vectors. The residual is defined by the signed distance between the points coming from the measurement systems after the change of reference frame and the stl mesh. This distance is calculated between the points of the closest scanning system and the STL mesh.
This problem being badly formulated, it is necessary to establish an initial solution in order to have a correct result.
This is why in a first step, we establish an initial using points that are identified on the STL model but also on the real model. A minimum of 3 points is necessary to block the 6DOF.
Then, in a second step, an optimization will be performed with the whole scan taking as initial solution the one previously calculated.
1.1.1.2. Modules importations¶
import trimesh
import numpy as np
from stl import mesh
from scipy import optimize
import plotly.graph_objects as go
import pandas as pd
import gbu
1.1.1.3. Case of study¶
Left Human scapula stl file
def stl2mesh3d(stl_mesh):
# stl_mesh is read by nympy-stl from a stl file; it is an array of faces/triangles (i.e. three 3d points)
# this function extracts the unique vertices and the lists I, J, K to define a Plotly mesh3d
p, q, r = stl_mesh.vectors.shape # (p, 3, 3)
# the array stl_mesh.vectors.reshape(p*q, r) can contain multiple copies of the same vertex;
# extract unique vertices from all mesh triangles
vertices, ixr = np.unique(
stl_mesh.vectors.reshape(p * q, r), return_inverse=True, axis=0
)
I = np.take(ixr, [3 * k for k in range(p)])
J = np.take(ixr, [3 * k + 1 for k in range(p)])
K = np.take(ixr, [3 * k + 2 for k in range(p)])
return vertices, I, J, K
def create_mesh3D(
stl_file="scapula.stl",
title="Mesh3d from a STL file<br>AT&T building",
colorscale=None,
):
if colorscale is None:
colorscale = [[0, "#e5dee5"], [1, "#e5dee5"]]
vertices, I, J, K = stl2mesh3d(mesh.Mesh.from_file(stl_file))
x, y, z = vertices.T
mesh3D = go.Mesh3d(
x=x,
y=y,
z=z,
i=I,
j=J,
k=K,
flatshading=True,
colorscale=colorscale,
intensity=z,
name="AT&T",
showscale=False,
)
layout = go.Layout(
paper_bgcolor="rgb(1,1,1)",
title_text=title,
title_x=0.5,
font_color="white",
width=800,
height=800,
scene_camera=dict(eye=dict(x=1.25, y=-1.25, z=1)),
scene_xaxis_visible=False,
scene_yaxis_visible=False,
scene_zaxis_visible=False,
)
fig = go.Figure(data=[mesh3D], layout=layout)
fig.data[0].update(
lighting=dict(
ambient=0.18,
diffuse=1,
fresnel=0.1,
specular=1,
roughness=0.1,
facenormalsepsilon=0,
)
)
fig.data[0].update(lightposition=dict(x=3000, y=3000, z=10000))
return fig
def add_points(points, fig, name, scale=1e3, *args, **kwargs):
fig.add_trace(
go.Scatter3d(
x=np.concatenate([points, points]).reshape(-1, 3)[:, 0] * scale,
y=np.concatenate([points, points]).reshape(-1, 3)[:, 1] * scale,
z=np.concatenate([points, points]).reshape(-1, 3)[:, 2] * scale,
name=name,
mode="markers",
**kwargs
)
)
return fig
fig = create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")
fig.show()
1.1.2. STL file¶
We define reference points on the STl that can be easily found on the real model.
A point on the glenoid, coracoid and acromion is selected.
p_glene_stl = np.array([0.02515289854664963, 0.03122873596800735, 0.036798809060995745])
p_coracoide_stl = np.array(
[0.026787610816186087, 0.05609095153197788, -0.0043762069034760384]
)
p_acromion_stl = np.array(
[-0.008420164259823573, 0.07140986464971243, 0.018059532550748165]
)
fig = add_points(points=p_glene_stl, fig=fig, name="glène")
fig = add_points(points=p_coracoide_stl, fig=fig, name="coracoide")
fig = add_points(points=p_acromion_stl, fig=fig, name="acromion")
fig.show()
The same points are easily found on the real model.
Here, the points are marked in green.
1.1.3. System scan device¶
1.1.3.1. Get scanned points from our system scan device¶
We get the points recorded by our system system scan device.
df_scan = pd.read_pickle("data_scan.p")
df_scan
| p_glene_scapula | p_acromion_scapula | p_coracoide_scapula | p_scan_scapula | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| p3d | p3d | p3d | p3d | |||||||||
| x | y | z | x | y | z | x | y | z | x | y | z | |
| 0 | 0.014046 | 0.010801 | 0.118422 | -0.023197 | -0.022039 | 0.086035 | -0.002085 | -0.037214 | 0.123242 | 0.015003 | -0.004692 | 0.11943 |
| 1 | 0.014152 | 0.01073 | 0.118276 | -0.023097 | -0.022086 | 0.085978 | -0.002091 | -0.037263 | 0.123189 | 0.015018 | -0.004673 | 0.119397 |
| 2 | 0.014185 | 0.010854 | 0.11838 | -0.023238 | -0.022144 | 0.086024 | -0.002113 | -0.037188 | 0.123366 | 0.015072 | -0.004562 | 0.119343 |
| 3 | 0.014362 | 0.01094 | 0.118554 | -0.023217 | -0.022149 | 0.086317 | -0.001946 | -0.037163 | 0.123109 | 0.015048 | -0.004463 | 0.119347 |
| 4 | 0.014379 | 0.010899 | 0.118494 | -0.023224 | -0.022023 | 0.086214 | -0.001959 | -0.036706 | 0.12319 | 0.015013 | -0.004507 | 0.11929 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 79 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.015566 | -0.001602 | 0.126646 |
| 80 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.01562 | -0.001514 | 0.126536 |
| 81 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.015541 | -0.001624 | 0.126567 |
| 82 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.015546 | -0.001497 | 0.126532 |
| 83 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0.015578 | -0.001497 | 0.126524 |
84 rows × 12 columns
1.1.3.2. Plot scanned points in their reference frame¶
Here the points are expressed in scapula reference frame
cols = set([cols[0] for cols in df_scan.columns])
fig = go.Figure()
for i, c in enumerate(cols):
fig.add_trace(
go.Scatter3d(
x=df_scan[c].p3d.x,
y=df_scan[c].p3d.y,
z=df_scan[c].p3d.z,
name=c,
marker_size=2,
mode="markers",
)
)
fig.show()
1.1.4. Create initial guess¶
1.1.4.1. Gathering input data of reference points¶
p_glene_scan = np.array(df_scan["p_glene_scapula"].p3d.mean())
p_coracoide_scan = np.array(df_scan["p_coracoide_scapula"].p3d.mean())
p_acromion_scan = np.array(df_scan["p_acromion_scapula"].p3d.mean())
tri_point_stl = np.concatenate([p_glene_stl, p_coracoide_stl, p_acromion_stl]).reshape(
-1, 3
)
tri_point_scan = np.concatenate(
[p_glene_scan, p_coracoide_scan, p_acromion_scan]
).reshape(-1, 3)
1.1.4.2. Compute initial guess¶
def unpack(X):
return X.reshape(2, 3)
def cost(X):
Rsc2stl, Tsc2stl = unpack(X)
tri_point_out = gbu.utils.apply_RT(tri_point_scan, Rsc2stl, Tsc2stl)
dist = abs(tri_point_stl - tri_point_out)
pen_pin_radius = 0
res = dist - pen_pin_radius
return res.flatten()
X0 = np.zeros(6)
sol = optimize.least_squares(
cost,
X0,
method="lm",
ftol=1.0e-12,
xtol=1.0e-12,
gtol=1.0e-10,
)
X_pre = sol.x
R_inital_guess, T_initial_guess = unpack(X_pre)
print(f"Initial guess solution: rvec={R_inital_guess}, tvec={T_initial_guess}")
Initial guess solution: rvec=[ 2.23440972 -0.14576349 0.55625112], tvec=[-0.02428643 0.12368749 0.09225085]
1.1.5. Surface match with full scan¶
1.1.5.1. Gathering input data of full scan¶
p_full_scan_sc = np.array(df_scan["p_scan_scapula"])
# load full scapula expand
scale = 1e-3 # convert to meter
mesh_scapula = trimesh.exchange.load.load("scapula.stl", type="stl")
mesh_scapula.vertices *= scale
1.1.5.2. Optimization¶
def cost(X, p3d, mesh):
Rsc2stl, Tsc2stl = unpack(X)
p_stl = gbu.utils.apply_RT(p3d, Rsc2stl, Tsc2stl)
dist = trimesh.proximity.signed_distance(mesh, p_stl)
pen_pin_radius = 0
res = dist - pen_pin_radius
return res
sol = optimize.least_squares(
cost,
X_pre,
method="lm",
ftol=1.0e-12,
xtol=1.0e-12,
gtol=1.0e-10,
args=(p_full_scan_sc, mesh_scapula),
)
R_sc2stl_sol, T_sc2stl_sol = unpack(sol.x)
print(f"Final solution: rvec={R_sc2stl_sol}, tvec={T_sc2stl_sol}")
Final solution: rvec=[ 2.2283109 -0.16742782 0.55011184], tvec=[-0.02341616 0.12482507 0.09038642]
1.1.5.3. Plot results¶
p_full_scan_stl = gbu.utils.apply_RT(p_full_scan_sc, R_sc2stl_sol, T_sc2stl_sol)
fig = create_mesh3D(stl_file="scapula.stl", title="Scapula Left <br> stl file")
fig = add_points(points=p_full_scan_stl, fig=fig, name="scanned_points", marker_size=10)
fig.show()
